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A novel scheme for the steady state solution of the standard Redfield quantum master equation is developed 
which yields agreement with the exact result for the corresponding reduced density matrix up to second 
order in the system-bath coupling strength. We achieve this objective by use of an analytic continuation 
of the off-diagonal matrix elements of the Redfield solution towards its diagonal limit. Notably, our scheme 
does not require the provision of yet higher order relaxation tensors. Testing this modified method for a 
heat bath consisting of a collection of harmonic oscillators we assess that the system relaxes towards its 
correct coupling-dependent, generalized quantum Gibbs state in second order. We numerically compare our 
formulation for a damped quantum harmonic system with the nonequilibrium Green's function formalism: 
we find good agreement at low temperatures for coupling strengths that are even larger than expected from 
the very regime of validity of the second-order Redfield quantum master equation. Yet another advantage 
of our method is that it markedly reduces the numerical complexity of the problem; thus allowing to study 
efficiently large-sized system Hilbert spaces. 



I. INTRODUCTION 

The application of canonical statistical mechanics in- 
herently assumes large environments interacting weakly 
with a few relevant degrees of freedom, then yielding the 
well-known canonical thermal state. The emergence of 
the thermal steady state canonical Gibbs density ma- 
trix at very weak coupling strength, or the generalized 
thermal Gibbs state at finite coupling strengths, when 
starting from quantum dynamical microscopic laws still 
presents a formidable problem. This objective is known 
under the label of open system quantum dynamics. A 
main goal is then to obtain the reduced system dynamics 
in terms of the reduced density matrix, which typically is 
approached using a wide variety of approximate quantum 
master equations. 

Formally exact generalized quantum master equations 
yield the reduced system dynamics either within a time- 
convolution (time-non-local) forroii^ and equivalently 
also in its time-convolutionless (time-local) form^Ti^. The 
hierarchy equations of motion approach^rJ^ yields yet 
another, formally exact approach in terms of an infi- 
nite number of auxiliary reduced density operators. All 
these formally exact approaches are computationally very 
demanding and, typically, can treat systems possessing 
a small Hilbert space dimension only. For these for- 
mally exact approaches it is only for specific setups, 
such as the situation involving (i) a system of harmonic 
oscillators^^"— , (ii) the intricate dissipative Landau- 
Zener dynamics at zero temperature^, or (iii) the known 
cases with a strictly pure dephasing dynamics^^r^^, that 
the exact solutions can be obtained. 



Timely applications, however, call for a definite need 
to study systems which span a rather large Hilbert space 
for its underlying nonlinear dynamics. In absence of 
analytic exact results the quantum master equations 
are typically evaluated using perturbation theory in the 
system(S)-bath(B) coupling strength. Commonly, the 
perturbation is truncated to second order in system- 
bath coupling, resulting in a whole group of approxi- 
mate quantum master equationa^^"— . Out of these many 
existing approximation schemes the Redfield quantum 
master equation (RQME) is the most generic one from 
which the Pauli^^ and the Lindblad^^i^ master equations 
can be deduced upon invoking further approximations^. 
Sometimes the RQME is also subjected to the secular 



approximatio: 
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and/or one neglects the Lamb shift- 
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type contributions^"—, which cannot always be justi- 
fied a priori. Generally, all these diverse approxima- 
tions, even within the weak coupling limit, do fail at zero 
temperature. This is so because of the neglect of alge- 
braic long-time tail contributions stemming from zero- 
temperature bath correlations^. Even without invoking 
such approximations, a question recently raised is the 
overall accuracy of the Redfield formalisnt^i^: Therein 
these authors demonstrated that the Redfield formalism 
is not correct for the steady state within its commonly 
used second order form. The discrepancy arises due to 
the second order diagonal elements which require contri- 
butions from the higher order relaxation tensor for their 
correct evaluation. In view of these findings it is not 
possible to capture correctly the effects of finite coupling 
up to second order by use of the Redfield formalism; - 
this feature also corrects some inadequately stated claims 
contained in the previous literature^i^. 

Our main goal with this study is to correctly evaluate 
the steady state reduced density matrix up to second or- 



der in the system-bath couphng without having to invoke 
higher order relaxation tensors. In order to successfuhy 
achieve this objective we put forward a modified solution 
of the Redfield quantum master equation by a procedure 
that uses the ofF-diagonal structure in second order to 
let approach its diagonal structure via a unique analytic 
continuation. 

Comparing this modified Redfield solution with rigor- 
ous canonical perturbation theor y "^°i"^^'^'^ , we show that 
the modified solution agrees with the exact reduced ther- 
mal equilibrium density operator; i.e. the generalized 
Gibbs state, reading: p = TrB(e~'^-f^'°')/Tr(e-'3^'°t), up 
to second order in the system-bath coupling strength. 
Our solution is not only accurate as compared to the 
RQME, but is also numerically efficient. This is because 
the inherent computational complexity in our method is 
of 0(A^'^), where N denotes the dimension of the sys- 
tem Hilbert space. Therefore, our technique enables 
one to quantum mechanically investigate the small-to- 
intermediate coupling strength regime for systems pos- 
sessing a large Hilbert space dimension. 

The paper is organized as follows: In Sec. |TT] we de- 
scribe our basic approach to model quantum dissipation 
and detail the RQME. In Sec. IIIII we elucidate the in- 
sufficient accuracy issue in the second-order steady state 
Redfield formalism. This part is followed by the expo- 
sition of our modified solution to the Redfield quantum 
master equation. In Sec. lIVI we consider a general nonlin- 
ear system that is connected to a harmonic bath and in 
the long-time limit show that it reaches the generalized 
Gibbs distribution within canonical perturbation theory 
carried out up to second order in system-bath coupling. 
In Sec. |V]we present the numerical comparison between 
our modified Redfield solution with the exact solution for 
a damped harmonic oscillator. We find a considerable im- 
provement between exact results and modified solution 
over extended regimes of weak-to-intermediate system- 
bath coupling strengths for which both the Redfield so- 
lution and the Lindblad solution fail. Sec. IVII summarizes 
our main findings while the Appendix details canonical 
perturbation theory. 



II. REDFIELD QUANTUM MASTER EQUATION IN 
PRESENCE OF ARCHETYPE QUANTUM DISSIPATION 

The basic approach to model quantum dissipation 
has been studied extensively before. The model Hamil- 
tonian for the bath and system-bath coupling has a 
long-standing history^i^"— but goes under the label of 
Zwanzig-Caledira-Leggett model^"— , 



Htat — Hs + Hg + i^RN + HsB, 



(1) 



denotes the generally nonlinear system Hamiltonian of a 
particle of mass M moving in a potential V{q). Here 
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describes the thermal environment as an infinite collec- 
tion of harmonic oscillators, each having a mass to„ and 
a frequency w„. 
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(4) 



is the potential renormalization in which the variable S 
denotes any function of the system variables p and q and 




^n-^n I 5 



(5) 



is the system-bath coupling Hamiltonian, wherein the c„ 
denotes the system-bath coupling constant of the n-th os- 
cillator with the system operator S. The collective bath 
operator is i? = — X]n=i '^nXn- Throughout this work we 
use ft = 1 and /cb = 1 . 

Using correlation-free initial conditions, i.e., Ptot(io) = 
Ps(io) (X" Puita), with Puita) being the canonical thermal 
state of the bath, and assuming overall weak system-bath 
coupling, we obtain the perturbative, 2-nd order Redfield 
quantum master equation in its time-local energy repre- 
sentation asSLSSi^i^, 
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(6) 



where the explicit time dependence in the reduced den- 
sity matrix p{t) = Tre (ptot(O) has been suppressed, i.e., 
Pnm = {'n\p{t)\my, |n) being the energy eigenvector of 
the bare system. Here, the matrix elements Sik are de- 
fined as Sik = \i|'5'|A:). Despite the apparent time-local 
form of Eq. (O the non-Markovian behavior is fully cap- 
tured due to the time dependence in the transition rates 
W. These are given by, 






jk 

7o 



W^, 



W',u 



where 
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Ej — Ek is the energy dif- 
v^J ference between system energy levels. Above 7o/2 which 



where VF'j, 



W'jk and A,fc 



arises from 7?^^ has been neglected from the uncoupled 
propagation (Eq. ([5|)) and accounted in W. 

Instead of specifying all the parameters of the bath we 
now define the spectral density J(w) as, 



3iu)^nJ2 
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Using the spectral density we can calculate the damping 
kernel at time i = 0; i.e., 70, used in Eq. ([7]) as. 
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and the equilibrium bath-bath correlation function 
C(t) = (^B{t)B^, where B{t) is evolving according to 
exp(— ii^BT), used in Eq. (|8]) as. 



C(r) = r ^J(.) 



coth ( — ) cos {lot) 



-i sin (wr) 



(11) 



Although the bath-bath correlation defined above is spe- 
cific to a harmonic bath, the theory presented here can 
readily be generalized to other bath models, e.g. spin 
baths as long as the bath-bath correlation C(t) can be 
evaluated. 



where A is a dimensionless parameter whose power in- 
dicates the corresponding order of the perturbation ex- 
pansion. Eventually, A will be set to 1. A above is a 
four tensor depending on the system Hamiltonian. The 
operator R'^")(t — i„) denotes the Redfield superoperator 
of rank 4 which depends both on the system operator 
and the bath correlators. Next we rearrange p into a col- 
umn vector and split it into its diagonal part (p^) and 
off-diagonal part (p^a)- Then, using the RQME (Eq. dH)) 
the 0-th order tensor in Eq. ([T2|) can be rewritten as a 
matrix assuming the form. 





A22 



(14) 



where A22 is a diagonal matrix with Aij (i ^ j) forming 
the diagonal. The four tensors R(")(i — t„) are also split 
accordingly; i.e.. 



R^"n^-io) 



i?il (i-to) R12 [i-ta) \ /-.f.^, 

4i'(i-to) 42^(t-to)r ^ ^ 



with no restrictions made for the form of the sub- 
matrices. For the specific case of n = 2, R(^)(t — t„) 
is the same as the Redfield tensor given in Eq. (j6l). 

In order to obtain the steady state we set dp/dt = 
and take the limits {t — t^) — > 00. Because the stationary 
problem is not dependent on time we will drop the paren- 
theses from the tensor, i.e., R^"^(oo) = R^"^. Therefore, 
using Eq. (fT2|) and Eq. (fTS]) we obtain: 



III. MODIFIED SOLUTION TO THE REDFIELD 
QUANTUM MASTER EQUATION 

A. Perturbative accuracy of steady state Redfield solution 

Recently Mori and Miyashita^ and, as well, Fleming 
and Cummings^ independently established that for a 
generic "2n" -order quantum master equation the solu- 
tion in the long-time limit can be correct only up to or- 
der "2n — 2" , because the diagonal elements loose their 
accuracy over evolving time. These authors suggest that 
in order to obtain a steady state solution correct up to 
2-nd order a 4-th order master equation^^^i^ should be 
used, which can be numerically accomplished for small 
system Hilbert spaces only. 

We first corroborate this finding with a different 
method; concentrating on the steady state accuracy of a 
2-nd order RQME. We start out with the generic pertur- 
bation series expansion to all orders in the system-bath 
coupling of the time-local, formally exact master equa- 
tion; i.e.. 



' 00 \ 00 

V n=2,4,6,--- / m=0,2.4,' 



A^p^") = 0. (16) 



In order to obtain p correct up to 2-nd order we equate 
the coefficients of the different powers of A to zero so that 
we obtain independent equations to calculate p*-*^' and 
pf^). This implies, 

1. Setting the co-efficient of A*^ equal to zero yields. 



pi? = 0. 



(17) 



2. Setting the co-efficient of A^ equal to zero implies. 
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3. Setting the co-efficient of A'^ equal to zero provides 
the condition. 



p(2) (2) _ _ p(2) (2) _ „(4) (0) 
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(20) 



^= A+ J2 ^"R^"Hi-io) 

\ n=2,4,6,-- > 

and the reduced density matrix, 

00 

p= J2 ^>^"^' 

n=0,2,4,--- 



(12) 



(13) 



Equation (|T9l) shows that in order to obtain the 2-nd 
order off-diagonal elements we need only the 0-th order 
and 2-nd order relaxation tensors which can be obtained 
from the RQME using Eq. ^. In contrast, in order 
to obtain the 2-nd order diagonal elements from Eq. ([20|) 
one requires knowledge of the 4-th order relaxation tensor 

-^11 ■ 



B. Analytic continuation procedure for diagonal density 
matrix elements 



In this section we present the procedure to obtain the 
stationary reduced density matrix that is correct up to 2- 
nd order in the system-bath couphng without the need to 
invoke the use of the 4-th order relaxation tensor. The 0- 
th order and the 2-nd order ofF-diagonal elements can be 
obtained correctly from the RQME as described above. 
Therefore, we use Eq. (fT7| and Eq. (fT8|) along with the 
Rcdficld tensor R(^) detailed in Eq. © to arrive at the 
0-th order reduced density matrix, 
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Therefore, Eq. ([23)) becomes, 
pI') = hm I - J2 SmS.n \(WU0) + K,i-^) 
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where. 



The 2-nd order ofF-diagonal elements follow from 
Eq. ((H]) as. 
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(n^m). (22) 



Note that if we next construct the diagonal elements by 
merely substituting n = m in Eq. (j22p . then the equation 
exhibits an indeterminate 0/0 singularity. This indicates 
that even though we cannot substitute n = m directly, 
the limit m ^ n might exist. If such a limit indeed exists 
and being unique, then by use of the uniqueness theorem 
the 2-nd order diagonal elements can be obtained by this 
limiting procedure. In order to perform this limit tti — >■ n 
we consider each element of the 2-nd order reduced den- 
sity matrix to be a function of the bare system energies 
Ei ii = I,--- ,iV). In the energy parameter space we 
vary only one of the energies Em and let it continuously 
approach the energy Em via a small complex parameter 
z; i.e., we set Em -^ En — z. 

In doing so, we start by splitting the transition rates 
in Eq. (|22p into its real and its imaginary parts, using 
Eq. ^ to obtain: 
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We next let £',„ — > En — z and perform the limit z — > 0. 
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Because prnm (being the un-normalized 0-th order re- 
duced density matrix) depends on the energy E^ we 

made use of the Taylor expansion of pmm around the 
energy En to retain up to the first order: 

UfJn 

We next define, 
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and note that \im,^oW^^{~z) = W/;(0) = W^^. 
Eq. p4|) can thus be recast as. 
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where, 



r7m 
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(w^„(o) + w;„(z))p; 



(29) 



In the limit z — > it follows from Eq. (P5|) that 

hm.^o W^;„;(-2)^ = lim,^oW^;;.(^) = W'nM = VF/„. 
Therefore, in this limit the term in the curly bracket in 
Eq. (pn|) assumes precisely the same form as Eq. (|2ip . 
hence it is equal to zero. Consequently, Eq. (|28p becomes. 
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Eq. ((30)) is independent of the way in which the energy 
Ejji approaches E^ and hence this hmit procedure is 
unique. The uniqueness of the hmit is crucial to ensure 
that the resuhing thermal steady state of the system is 
unique. 

The diagonal elements of the density matrix obey the 
normalization condition Tr(p) = 1. Since we performed 
an analytic continuation to obtain the 2-nd order diag- 
onal elements there is no guarantee the normalization 
condition is preserved. Therefore we can write the nor- 
malization condition explicitly as, 
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where we have ignored the 4-th and higher order terms 
and used the condition ^^ pi^ = 1 , which is required 
to determine p^^' uniquely. Therefore, upon normalizing 
Eq. (pO|) with help of Eq. ([3T|) we obtain the first main 
result, 
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Using Eq. ([5^ to calculate the 2-nd order diagonal ele- 
ments we need to know the derivative of the 0-th order 
reduced density matrix, dpnn/dEn- This derivative de- 
rives from Eq. (HI]), which is satisfied by p^°^ and sub- 
sequently differentiate with respect to the energy En to 
find 



dp)iA 2^i^n^rii^i''i\-dA~Pii ^ ^K~ P 
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Therefore, we have all the ingredients at hand to cal- 
culate the 2-nd order diagonal elements from Eq. ([5^ : 
This constitutes the first main result of our work. In 
our derivation we have made no assumptions besides the 
validity of analytic continuation. The above outlined the- 
ory can be readily generalized to multiple heat baths; a 
topic to be addressed by us in future work^. 

The modified solution outlined above is correct not 
only up to 2-nd order in system-bath coupling but addi- 
tionally it is well suited for numerical studies: Numerical 
simulations with the RQME are very cumbersome be- 
cause the relaxation tensor R'^^ in Eq. ((6)) scales as the 
fourth power— of the system Hilbert space dimension N. 
Therefore, in the steady state the computational com- 
plexity of the problem typically scales proportional to 
N^, assuming that the analytic forms of the transition 
rates are known. On the other hand, in our modified 
Redficld solution all components of the reduced density 
matrix can be obtained by reference to the transition 



rates W only, which scale as 7V^. Thus, in the modified 
solution the computational complexity becomes drasti- 
cally reduced to be of order N^ . This fact is equivalent 
to solving the quantum master equation with use of the 
continued fraction scheme^; it thus enables us to study 
systems with much larger Hilbert space dimension. 



IV. COMPARING MODIFIED REDFIELD SOLUTION 
WITH SECOND ORDER CANONICAL PERTURBATION 
THEORY 

For a finite system-bath coupling the thermal equilib- 
rium density matrix is typically no longer of Gibbs type 
(strict weak coupling limit) but rather of the general- 
ized Gibbs form, p°'^ ex TrB(e~''-^'°'), resulting in a quan- 
tum Hamiltonian of mean forced. It is interesting to 
know if this distribution can be obtained from a full non- 
Markovian dynamical theory of a system weakly coupled 
to a heat bath. Although this seems reasonable there 
is no agreed consensus on this issue from the viewpoint 
that the literature deals with a variety of perturbative 
quantum (2-nd order) master equations^^il^Si^. 

Since the Redficld formalism is rigorously valid only in 
the A — > limit it is expected that in this very limit the 
canonical form p°^ (x e~^^^ emerges. In order to test 
the accuracy of our novel modified Redficld solution we 
implement an order by order comparison between canon- 
ical perturbation theory (CPT), which perturbatively ex- 
pands the generalized Gibbs distribution, as detailed in 
the Appendix, with our modified Redficld solution. Ac- 
cording to CPT (Eq. (|M3)) . dHH), and (|M5|) ) the re- 
duced density matrix up to 2-nd order in the system-bath 
coupling reads 
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wherein the different contributions assume the form: 
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A. Comparing the 0-th order result 

Let us first compare the 0-th order reduced density 
matrix. For the harmonic baths described by Eq. (jH]) it 
can be shown that the bath-bath correlator C(r) obeys 
the Kubo-Martin-Schwinger (KMS) conditioii^^i^, 



C{-T)=C{T-il3). 



(38) 



This implies that the real part of the transition rates W' 
obey the detailed balance condition^ given by, 
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The analytic form of the 0-th order reduced density ma- 
trix can be obtained upon using Eq. (|2ip as, 
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where Z^ = ^i e~" '- . A direct comparison between 
Eq. (|40|) and Eq. p4|) yields the expected result that at 
the 0-th order CPT agrees with our modified, 0-th order 
Redfield solution. 



B. Comparing the 2-nd order result 

More intriguing is the comparison of the modified Red- 
field solution with the 2-nd order CPT-result. The 2-nd 
order reduced density matrix obtained from CPT can be 
manipulated further so that it indeed matches precisely 
our modified Redfield solution. In order to demonstrate 
this fact wc first simplify the integral occurring in £), 
Eq. (|36l). by using the definition of the bath-bath corre- 
lator C(r) in Ec^. ([TT|) to obtain 
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(41) 

where we have interchanged the w- (stemming from C(t)) 
and X- integration and performed the x-integral analyt- 
ically. We next express the right hand side in terms of 
the transition rates W, which enter in our modified Red- 
field solution. In order to do this we use the so termed 
Sokhotskyi-Plcmelj formula^, 
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Here, P denotes the principal value. Therefore, using the 
above identity along with Eq. ([8]) and Eq. ([TT|) we can 
express the imaginary part of the transition rates W" in 
the form 
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Therefore using the above equation, Eq. (|4T|) can be ex- 
pressed as. 



da;e-^^'^ C{-ix) = W^'.+c-^^^^ W'',. (44) 



1. OfT-diagonal elements 

Upon use of Eq. (|44)) the 2-nd order off-diagonal ele- 
ments from CPT, i.e., Eq. ([SS]) can be expressed in terms 
of W" as 
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where we have absorbed the 70 into W" , according to 
Eq. (I7|) . Formally adding the real part of the transition 
rates W' into Eq. (|^5|) , but noting that this so added con- 
tributions vanish identically by virtue of detailed balance 
in Eq. ([55]) . we find the result 
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Upon comparing Ec^. ^F^ with Eq. (|46p wc find that the 
CPT and our modified Redfield solution arc identical. 



2. Diagonal elements 

Most importantly, we next test the agreement between 
the 2-nd order diagonal elements from CPT with our 
modified Redfield solution. Noting that the integral oc- 
curring in Eq. ([55]) is the derivative of Eq. pi)) w.r.t A^ 
we obtain, 
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where Vij has been defined in Eq. (P7|) . Because 

dp^°'^/dE, = -I3p^°\ Eq. dSni) exactly matches Eq. gH). 
This constitutes a second main result: Namely CPT up 
to 2-nd order and our modified Redfield solution are in- 
deed in perfect agreement. This shows that in the weak, 
but finite coupling limit the long-time thermal reduced 
density matrix stemming from a non-Markovian theory 
is of the generalized Gibbs form. 



V. A TEST CASE: 
OSCILLATOR 



DAMPED HARMONIC QUANTUM 



In this section we compare our modified Redfield 
solution to the exact nonequilibrium Green's function 
(NEGF) resulta^S for a damped harmonic oscillator that 
is linearly coupled to a thermal heat bath. This compar- 
ison will allow us to estimate the system-bath coupling 
strengths that can be probed safely by employing the 
presented modified scenario. In order to do so we use 
a specific spectral density 3{uj) of the heat bath. Sev- 
eral phenomenological forms of the spectral density are 
in use in the literature for the damped oscillator quan- 
tum dynamics; sometimes a numerical decomposition is 
employed to save computational costs^i^. Here we use 
the common Lorentz-Drudc (LD) spectral density, 



J(w) 



M^LO 



1 + {uj/lu-d) 



(48) 



where u-o denotes the cutoff frequency and 7 is the phe- 
nomenological Stokesian damping coefficient which char- 
acterizes the system-bath coupling strength. Since 7 oc 
S^i c^i the spectral density is of 2-nd order in system- 
bath coupling. Decomposing the hyperbolic cotangent in 
Eq. (fTTj) into its Matsubara frequencies, vi = 2TrlT, where 
T denotes the temperature of the bath, and noting that 
the resultant equation exhibits poles at a; = ±i Wd and 
uj ~ ^ii^i we can calculate C(t) explicitly by use of the 
residue theorem to obtain 



C(t) = ^w^cot 



/3a;o 



2Mj 



St 



Vie 



1=1 



. M-f 



ujI c-"°^ 



sgn(r). 



{vi/uj^Y 
(49) 



Therefore, the components of the transition rates W , de- 
fined in Eq. ([5]), read 



W 
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(50) 
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with the relation that 



7o = 7Wd- 



(51) 



(52) 



The system Hamiltonian is a single harmonic oscillator, 
reading 



H, 



^^l"-'-'- 



(53) 




FIG. 1. (Color Online) Plot of the discrepancy error DE'^ of 
the ground state population versus the dimensionless system- 
bath coupling strength (7/1^0) for a damped harmonic quan- 
tum oscillator. Top panel: The (black) solid line depicts the 
rather small discrepancy for our modified Redfield solution 
{X = MRS) and the bottom panel (red) dashed line shows 
the large discrepancy obtained via the ordinary Redfield quan- 
tum master equation {X = RQME). Our parameters used for 
the calculation are M = 1 u, ljo = 1.3 x lO" Hz, T = 50 K 
and the cutoff is chosen at ljd ~ 10 ojq. 



where x, p, M, and ujq are the position, momentum, 
mass and angular frequency of the oscillator. The har- 
monic oscillator is linearly coupled to the bath via the x- 
coordinate. This implies that S = x in Eq. ([5]). Through- 
out this work the system-bath coupling will be measured 
in a dimensionless parameter, defined by taking the ratio 
between the damping coefficient 7 and the angular oscil- 
lator frequency, i.e., 7/wo- Since the Redfield formalism 
is formulated in terms of cigcnbasis of the system Hamil- 
tonian of finite dimension, we choose a system Hilbert 
space that is sufficiently large so that even at the highest 
temperatures the occupation probability of finding the 
particle in the highest available energy levels is practi- 
cally zero. We do this by iteratively increasing the size 
of the system Hilbert space until at least five largest en- 
ergy levels possess a population less than 10~^^: In our 
case of the damped harmonic oscillator this results in 
around 40 levels. Using these 40 levels we can cover a 
temperature range up to five times the Debye tempera- 
ture, Td = {huJo)/kB. 

The main goal in this work is to correctly evaluate the 
2-nd order diagonal elements. At the 0-th order level 
the RQME and our modified Redfield solution give the 
canonical solution, which matches the result obtained 
from the NEGF method by taking the zero coupling limit. 
Therefore, in order to sensitively compare the 2-nd order 
elements we define a relative discrepancy error DE as 
follows: 



DE^ 



h/^o) 



(54) 



where p^^'^^ denotes the exact reduced density matrix 
obtained from NEGF method, p^ is the reduced density 
matrix obtained from the perturbative method; being ei- 
ther our modified Redfield solution (X = MRS) or the 
Rcdficld quantum master equation (X = RQME), and 
the ratio 7/0^0 specifies the overall system-bath coupling 
strength. Since the 4-th order term of the reduced den- 
sity matrix is of the order of {^/uJq)'^, it is expected that 
the discrepancy error is one order lower, i.e., 0{-f/uJo) 
if and only if the 2-nd order elements are calculated cor- 
rectly. In order to check this behavior we plot the discrep- 
ancy error versus 7/1^0 for the ground state population 
at T = 50 K in Fig. [H 

Since the temperature is chosen low the population 
Pii presents an appropriate quantifier for the complete 
reduced density matrix. The figure depicts that the dis- 
crepancy error for our modified Redfield solution (solid 
black line) indeed stays throughout of the order of (7/wo) 
for all coupling strengths 7/cJo; in contrast, for the 
RQME (dashed red line) the discrepancy error grows in 
absolute value 3> O (7/(^0 ), indicating the inaccuracy in 
the 2-nd order elements. In the limit {j/loo — ^ 0) the dis- 
crepancy error should vanish: This holds true only for our 
modified Redfield solution whereas the Redfield solution 
depicts a finite value which indicates that this 2-nd order 
(i.e., being proportional c^) solution indeed is not correct 
to leading 2-nd order. The temperature does not play a 
major role, since all the features in the discrepancy error 
remain the same up to temperatures of ^ 3000 K. There- 
fore, for all values of the weak-to-moderate system-bath 
coupling strengths our modified Redfield solution is able 
to predict the 2-nd order elements correctly, whereas the 
RQME fails to do so already for small values of (7/^0)- 

In Fig. [2] we study the actual population values for the 
first few levels as a function of 7/ajo for two different tem- 
perature. Fig. Wia) corresponds to a temperature of T = 
50 K and[2i;b) is for T = 1000 K. We have opted to plot 
these two extreme temperatures because for all interme- 
diate temperatures values the features of the plot remain 
practically the same. We compare our modified Red- 
field solution (black solid line) with the RQME (dashed 
red line), NEGF results (blue crosses) and the Lindblad 
master equation (dotted green line). The Lindblad mas- 
ter equation is extensively used in the litcratur o^^'^^ due 
to its ease in computation and its preservation of pos- 
itivity. Although positivity is an essential criteria, the 
Lindblad solution for the damped quantum harmonic os- 
cillator case is the canonical distribution^ with no ex- 
plicit dependence on coupling strength. Put differently, 
the Lindblad solution always fails to capture the effects 
of finite system-bath coupling. On the other hand the 
RQME depicts severe deviations from the exact result 
for small, but finite coupling strengths. 

In the extreme low temperature regime the RQME 
yields negative populations, note the results for (/CI22, P44) 
in panel Fig. HJa) already for weak coupling strengths, 
indicating that the validity of the solution holds only in 
the zero coupling limit. The steady state solution of the 
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FIG. 2. (Color Online) Graph of the populations for the 
first four lowest lying energy levels versus the dimensionless 
system-bath coupling strength 7/0)0 for a damped quantum 
harmonic oscillator. The (black) solid lines correspond to our 
modified Redfield solution, the (red) dashed lines present the 
results for the Redfield quantum master equation, the (green) 
dotted-dashed lines depict the results for the Lindblad solu- 
tion, while the (blue) crosses represent the exact result using 
NEGF. Panel (a) is for the temperature of T = 50 K and 
panel (b) corresponds to a temperature of T = 1000 K. The 
remaining parameters used for the calculation are M = 1 u , 
ojo = 1.3 X 10^'' Hz and cjd = 10 oja- 



RQME has been critiqued befor o^^'^^ for producing un- 
physical, negative populations. We can now assess that 
the reason for its breakdown is rooted in the incorrect 
2-nd order diagonal elements. The modified Redfield so- 
lution matches the exact solution quite well for system- 
bath couplings 7/ci;o as strong as 0.2, even at low tem- 
peratures. At high temperatures our modified solution 
yields a most impressive agrement with the exact results, 
extending over sizable regimes of coupling strengths up 
to ^/ujg > 0.6. Beyond a coupling strength 7/0;,, ~ 0.6 
the modified Redfield density matrix is no longer pos- 
itive definite, which is determined upon examining the 
eigenvalues of the reduced density matrix; this indicates 



a breakdown of 2-nd order perturbation theory beyond 
this value. Nevertheless, the presented modified Redfield 
solution provides a decisive and salient improvement over 
the RQME in that the coupling strengths that can be 
probed accurately becomes sizable. 



VI. CONCLUDING REMARKS AND OUTLOOK 

In this paper we have demonstrated via the exact com- 
parison with canonical perturbation theory and extensive 
numerics that the Redfield quantum master equation is 
inaccurate for the steady state. This failure is the re- 
sult of incorrect second order diagonal elements. This is 
in the spirit of remarks made before by Fleming et al^. 
Their suggestion in overcoming this flaw by the use of the 
fourth order tensor to improve the second order accuracy 
is numerically extremely cumbersome. Their suggestion 
to use instead the Davies approximation^^— in order 
to obtain the thermal steady state reduced density ma- 
trix correct up to second order is not appropriate either. 
This is so because in case of the Davies approximation 
the long-time limit can be taken only if we take A — )• 0, so 
that the product \^t remains constant^. This in turn im- 
mediately implies that the Davies approximation is pre- 
cise only to 0-th order accuracy in the long-time limit, 
where it agrees with the Lindblad solution. Attempts 
have been made to correctly evaluate the second order 
diagonal elements using the Dyson expansion^S: In this 
context it must be noted, however, because the Dyson 
series is asymptotically divergent^, and although the off- 
diagonals in fact agree, the second order diagonal ele- 
ments are found not to match the exact result of Dhar 
et al^ for the damped harmonic oscillator problem. 

Therefore, since most of the perturbative methods fail 
to capture the effects of finite system-bath coupling in the 
long-time limit, we put forward a modified solution to the 
Redfield quantum master equation which reproduces the 
second order elements exactly. The derivation is based 
on obtaining the second order diagonal elements from 
the off-diagonal ones using an analytic continuation pro- 
cedure as detailed in Sec. IIIIBI The result of this scheme 
is unique indicating a unique steady state for the reduced 
density matrix. In order to test the validity of our solu- 
tion, we have compared our modified Redfield solution 
to canonical perturbation theory and demonstrated that 
our modified solution agrees with the generalized Gibbs 
distribution up to second order in the system-bath cou- 
pling strength for a general system that is coupled to a 
harmonic oscillator bath. This indicates that even in the 
weak, but finite, system-bath coupling limit the system 
thermalizes to a generalized Gibbs distribution. As will 
be elaborated elsewhere our method is also applicable to 
systems connected with multiple-baths, thus exhibiting 
nonequilibrium steady state transport^. 

As an illustrative example we tested and compared in 
Sec. |V] the reduced density matrix obtained by our mod- 
ified solution, the Redfield formalism, and the Lindblad 



master equation against the exact NEGF results for a 
damped harmonic oscillator. We find that our modified 
solution agrees quite well with the exact result for cou- 
pling strengths as strong as 7/wo = 0.2, showing a major 
improvement over the RQME which matches the exact 
result only in the limit j/luo — > 0. On the other hand the 
Lindblad solution for the damped oscillator case always 
yields the canonical distribution, wrongly indicating that 
the solution is not affected by the system-bath coupling 
strength. 

The presented modified Redfield solution further is nu- 
merically very efficient; this is mainly so because with 
our scenario the computational complexity scales as N^, 
where A^ is the system Hilbert space dimension, as com- 
pared to N^ for the Redfield formalism. This fact allows 
us to describe accurately not only the effects of finite 
system-bath coupling, but as well as to explore systems 
with rather large Hilbert space dimensions. A yet un- 
solved challenge consists in the extension of our scheme to 
the time- dependent relaxation of the reduced density ma- 
trix p(t) and; in this context, also the extension to study 
the differing relaxation processes that stem from different 
initial preparation schemes away from the typically used 
case of a correlation-free initial preparation. Assuming 
bath spectral densities that assure an ergodic behavior, 
the long-time limit is not affected by the initial prepara- 
tion, being in distinct contrast to its temporal relaxation. 
Yet another unsolved objective presents the perturbative, 
accurate study of multi-time correlations of open system 
observables, both time-homogeneous thermal and time- 
dependent nonequilibrium correlations beyond the weak 
coupling limit. 
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APPENDIX: CANONICAL PERTURBATION THEORY 

With this Appendix we outline the basic reasoning 
underlying canonical perturbation theory^i^i^ (GPT). 
This will assist us in determining the correct equilibrium 
reduced density matrix up to second order in coupling 
strength for a harmonic bath which is coupled bi-linearly 
to a general system. The basic idea dates back to the 
works of Peierls^ and Landau^ who calculated the free 
energy of the full system using a similar expansion. Here 
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we employ similar techniques for the reduced density ma- 
trix, which in case of the equilibrium problem is well de- 
fined by the generalized Gibbs distribution^: 



P 



cq 



TrB c-'3^t= 



(Al) 



where iJtot is defined in Eq. ([T]). We now use the Kubo 
identity^, 



p/3(A+S) _ pPA 



dAe'^^Be^(^+^) 



(A2) 



which is exact. Upon expanding e"*^^'"' up to second 
order in the coupling strength and taking the trace over 
the bath degrees of freedom we obtain 



TrB(c-^^'°') 



-l3Hs 



7o 



d/3i^(-i/3i)5(-i/3i) 



+ d/3i / d/325(-z/3i)5(-z/32)C(-i(/3i-/?2)) 
Jo Jo 

(A3) 

where S{—i(3i) = e'^i-^s g q-PiHs jg ^i^^ -fj.gg evolving sys- 
tem operator in imaginary time and C(— i(/3i — /32)) is 
the imaginary-time bath correlator as defined in Sec. |lll 
Using Eq. ^Ml in Eq. ^M^ the CPT reduced density 
matrix thus reads 



-PHs £, e-^^^TrsCl?) 



(ZsY 



(A4) 



where, 



D= d/3i / d/32^(-^/3l)^(-^/32)C(-i(/3l-/32)) 
Jo Jo 



7o 



2 Jo 



d/3i5(-i/3i)5(™i/3i). 



(A5) 
(A6) 



Next writing Eq. (jA4|) in the basis of the system Hamil- 
tonian we obtain. 
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wherein 



d/SiC^i'^"' / 'd/32e^^^-"C(-i(/?i-/32)) 
Jo 



7o 



A* 



d/3ie' 



,;3iA„„ 
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In Eq. dMI) A, 
in Eq. (gl). 

The main task in CPT is to evaluate the elements of 
the matrix D, Eq. (jASp . In order to do this we split the 
matrix D into its diagonal and off-diagonal elements and 
deal with each part separately, as detailed below. 



A. OfF-diagonal elements of the matrix D^m 

In order to obtain the off-diagonal elements of the 
matrix D we make the following change of variables: 
X = Pi — 132, y = Pi + 1^2 and then we perform the y 
integral analytically to find 



D^ 



— Y, (DnlSlm - DrnlSl„) , (A9) 



where, 



Dni = Sr^i e-P^- y dxC(-z x) e--^- -| 



B. Diagonal elements of the matrix !?„„ 



(AlO) 



For the diagonal elements of Z?, by using the same set 
of transformations as before, the integrals simplify and 
the diagonal elements of matrix D emerge as 



D^ 



where, 

Dnl = Snl C-^^" 
/9 



; 
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dxC{—ix)c 



-xA,, 



To 
2 



da;C(-ix)a;e-^^'" 



(A12) 



In summary, the thermal equilibrium reduced density 
matrix obtain via CPT is given, up to second order, by 
the generalized Gibbs state, reading: 

(A13) 



„CPT ^ (0),CPT , (2), CPT 

rnni rnm ' rnm 



where. 



(0),CPT 

rnm 



(2),CPT 
rnm 



-I3E„ 



D„ 



,-/3-E„ 



E.A. 



z. 



{Zsf 



(A14) 
(A15) 



Here, the off-diagonal elements of Dnm are given by 
Eq. (|A9p and the diagonal elements are given by 
Eq. (|XTT|) . ((IT2|) . Eq. ([MS]) exhibits that the equilib- 
rium reduced density matrix obtained via CPT is hcr- 
mitian and is normalized properly with trace over the 
system degrees of freedom equal to 1 . 
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